function M = mass_bb_2(dofs, V, T, d, varargin)
%        M = mass_bb_2(dofs, V, T, d, varargin)
% This is hp version, where d is a vector
%

if nargin == 4    
    d_max = max(d);
    Mat = cell(d_max,1);
    for degree = 1:d_max
        Mat{degree} = build_dim2(degree,degree);
    end
else
    Mat = varargin{1};
end

%% calculate areas for each triangle
x21 = V(T(:,2),1)-V(T(:,1),1);
y21 = V(T(:,2),2)-V(T(:,1),2);
x31 = V(T(:,3),1)-V(T(:,1),1);
y31 = V(T(:,3),2)-V(T(:,1),2);
areas = 0.5*(x21.*y31 - x31.*y21);

m = 0; nt = size(T,1);
for k=1:nt
    degree = d(k);
    md = (degree+1)*(degree+2)/2;
    m = m + md^2;
end
Indx1 = zeros(m,1);
Indx2 = zeros(m,1);
S = zeros(m,1);
pos_start = 1;
% tic;
for tri = 1:nt
    degree = d(tri);
    [i,j,s] = find(areas(tri)*Mat{degree});
    loc_dof = dofs(tri,1): dofs(tri,2);
    L = length(i);
    Indx1(pos_start:(pos_start + L-1)) = loc_dof(i);
    Indx2(pos_start:(pos_start + L-1)) = loc_dof(j);
    S(pos_start:(pos_start + L-1)) = s;
    pos_start = pos_start + L;    
end
% toc;
dim = max(max(dofs));
M = sparse(Indx1(1:(pos_start-1)),Indx2(1:(pos_start-1)),S(1:(pos_start-1)),dim,dim);
